#include <iostream>
#include <cmath>
#include "../sparse_matrix/sparse_crs.h"

const double PI = (std::atan(1.0) * 4.0);

double source(double x, double y);

double u_b(double x, double y);

double source(double x, double y)
{
    return 2.0 * PI * PI * std::sin(PI * x) * std::sin(PI * y);
};

double u_b(double x, double y)
{
    return std::sin(PI * x) * std::sin(PI * y);
};

int main(int argc, char *argv[])
{
    double x0 = 0.0;
    double y0 = 0.0;
    double x1 = 1.0;
    double y1 = 1.0;

    int N = 1000;
    double h = (x1 - x0) / (N - 1);
    int m = N - 2;
    int n = m * m;

    struct PATTERN pA;	/**< 系数矩阵模板 */
    struct MATRIX A;	/**< 系数矩阵 */
    
//    double f[n];
    double *f = new double [n];

    init_pattern(&pA, n, n, 5);

    for (int i = 0; i < n; i++)
	f[i] = 0;

    for (int i = 0; i < m; i++)
	for (int j = 0; j < m; j++)
	{
	    int idx = i * m + j;
	    set_pattern(&pA, idx, idx);
	    if (j != 0)
		set_pattern(&pA, idx, idx - 1);
	    if (j != m - 1)
		set_pattern(&pA, idx, idx + 1);
	    if (i != 0)
		set_pattern(&pA, idx, idx - m);
	    if (i != m - 1)
		set_pattern(&pA, idx, idx + m);
	}
    compress_pattern(&pA);
    init_matrix(&A, &pA);
    for (int i = 0; i < m; i++)
	for (int j = 0; j < m; j++)
	{
	    int idx = i * m + j;
	    double lx = (j + 1) * h;
	    double ly = (i + 1) * h;
	    add_matrix(&A, idx, idx, 4.0);
	    f[idx] = source(lx, ly) * h * h;
	    if (j != 0)
		add_matrix(&A, idx, idx - 1, -1.0);
	    else
		f[idx] += u_b(lx - h, ly);
	    if (j != m - 1)
		add_matrix(&A, idx, idx + 1, -1.0);
	    else
		f[idx] += u_b(lx + h, ly);
	    if (i != 0)
		add_matrix(&A, idx, idx - m, -1.0);
	    else
		f[idx] += u_b(lx, ly - h);
	    if (i != m - 1)
		add_matrix(&A, idx, idx + m, -1.0);
	    else
		f[idx] += u_b(lx, ly + h);
	}
//    double u[n];
    double *u = new double [n];
    for (int i = 0; i < n; i++)
	u[i] = 1.0;
    double err = 1.0;
    int k = 0;
    double tol = 1e-12;
    int maxit = 100000;
    while (err > tol)
    {
	gs_step(&A, u, f);
	err = res_norm(&A, u, f) / n;
	k++;
	if (k > maxit)
	{
	    printf ("not convergent!\n");
	    break;
	}
    }
    

    std::cout << "u = [" << std::endl;
    for (int i = 0; i < n; i++)
	std::cout << u[i] << std::endl;
    std::cout << "];" << std::endl;
    delete [] f;
    delete [] u;
    return 0;
};
